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A geometric analysis of the SHAKE and RATTLE methods for constrained 
Hamiltonian problems is carried out. The study reveals the underlying dif- 
ferential geometric foundation of the two methods, and the exact relation be- 
tween them. In addition, the geometric insight naturally generalises SHAKE 
and RATTLE to allow for a strictly larger class of constrained Hamiltonian 
systems than in the classical setting. 

In order for SHAKE and RATTLE to be well defined, two basic assumptions 
are needed. First, a non-degeneracy assumption, which is a condition on 
the Hamiltonian, i.e., on the dynamics of the system. Second, a coisotropy 
assumption, which is a condition on the geometry of the constrained phase 
space. Non-trivial examples of systems fulfilling, and failing to fulfill, these 
assumptions are given. 
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1 Introduction 

SHAKE and RATTLE are two commonly used numerical integration methods for Hamilto- 
nian problems subject to holonomic constraints [16, 1, 11, 9, 15]. The difference between 
the two methods is that RATTLE preserves "hidden" constraints, whereas SHAKE does 
not. For details and a historical account, see the monographs [10, §7.2] and [8, § VII. 1.4]. 

In this paper we give a rigorous geometric analysis of the SHAKE and RATTLE meth- 
ods. Our approach is based on an observation by Reich [15], that SHAKE and RATTLE 
may be expressed using flow maps. The analysis sheds light on the underlying "geo- 
metric foundation" of the two methods, and also on the exact relation between them. 
In addition, the geometric insight allows us to integrate a larger class of constrained 
problems than before. Indeed, the geometric versions of SHAKE and RATTLE work for 
coisotropic constraints. This class of constraints is strictly larger than the class of holo- 
nomic constraints. In particular, they may depend on both position and momentum. 

Throughout the paper we utilise the language of differential geometry. The main 
reason for doing so is not to generalise SHAKE and RATTLE from R 2d to manifolds, but 
rather because this notation makes the geometric structures more transparent. However, 
in order to link to the standard literature on SHAKE and RATTLE, we give many key 
results also in the classical K? d setting as examples. 

Our notation closely follows that of Marsden and Ratiu [13]. In particular, if Ai 
and J\f are two manifolds and / S C°° (M. , M) , then T/ : TM. — > TJ\f denotes the tangent 
derivative. If M C M. is a submanifold, then ij\f. M — > M. denotes the inclusion, and 
the pullback of differential forms from M to M is denoted i\f. If (V,(jo) is a symplectic 
manifold, and H 6 C 0O (7 :> ), then Xh denotes the corresponding Hamiltonian vector field, 
and the Poisson bracket is denoted {•,•}. The contraction between a vector X and a 
differential form a is denoted \xa. 

We continue this section with an outline of the paper and the main results. 



2 



Problem formulation Let (V,uj) be a symplectic manifold, H G C°°{V) a smooth 
function, and M C P a submanifold. Given (V,uj,H,M), the problem is to find a 
smooth curve 1 1- )• 7(f) such that 

^(Lyw-dtf)=0, 7(t)eX. (la) 

Equation (la) looks like a Hamiltonian system on but constrained to stay on the 
submanifold A4, called the constraint manifold. It can be rewritten as 

i^u - dH M = 0, (lb) 

where v = i^oj and Hm = 1 m^- From this formulation it is clear that equation (1) is 
intrinsic to A4, i.e., it only depends on objects defined on A4. 

Example 1.1. Let V = H 2d with canonical coordinates z = (q,p), and let A4 = 
ff _1 ({0}), where g = (gi,...,g m ) G C°°(R 2d ,R m ). Equation (1) then takes the form 

q = H p (q,p) +g p (q,p) T X 

P = -H q (q,p) - g q (q,p) T X (2) 
= g{q,p), 

where g q and g p are the partial Jacobian matrices, and A = (Ai, . . . , A m ) are Lagrange 
multipliers. Notice that if g does not depend on p, then this reduces to a canonical 
Hamiltonian system with holonomic constraints. 

Owing to the Hamiltonian structure of the equations, the reduction of the differential- 
algebraic equation (DAE) (2) to an ordinary differential equation (ODE) takes a particular 
form. This was already noticed by Dirac [2, § 1], and later perfected in [4]. 

Existence and uniqueness Since equation (1) is intrinsic to Ai, it is clear that any 
condition or assumption for existence and uniqueness should also be intrinsic, so it is 
enough to investigate existence and uniqueness intrisically on A4. 

The 2-form v is closed, but in general degenerate, so is not, in general, a 

symplectic manifold. Instead, it is a presymplectic manifold. 

The kernel of v defines a distribution on M denoted kerz;. As detailed in § 2.1, the 
kernel distribution is integrable. Thus, for each z G M, there is a submanifold K. z C M. 
such that T x /C 2 = ker v x for each x G fC z . 

If 7(t) is a solution to (lb), then dHj^^^)) G w(T '^^M). Thus, solutions stay on 
the hidden constraint set, given by 

M H = {z G M : dH M (z) G v(J z M)}. (3) 

In particular, a necessary condition for equation (lb) to have a solution is that the initial 
data fulfils 7(0) G M H , which is assumed from here on. 

As is further explained in §2, a sufficient condition for (local) existence and uniqueness 
of equation (lb), and hence equation (1), is the following. 
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Assumption 1 (Nondegeneracy) . For any z £ Ai, the critical points of the function 
Hjc x = i*fc^H_M are nondegenerate. 

Example 1.2. For the classical setting in Example 1.1, 

M H = {z £ R M : 9i (z) = 0, X H (z) ■ V 9i (z) = 0, * = 1, . . . , m} 

and Assumption 1 means that the matrix g z (z) T H zz (z)g z (z) is invertible for z £ M H . 
If g does not depend on p, then this is slightly weaker than the classical assumption that 
9qHppg g is invertible (see Remark 5.1). 

Geometric SHAKE and RATTLE We now define SHAKE and RATTLE geometrically. 

Definition 1.3 (Geometric SHAKE). Let cph be a method approximating exp(hXn)- 
The SHAKE mapping : Ai 3 zq \-t Z\ £ Ai is defined implicitly by 

<Ph(y) 6 M, y G 1C Z0 n O, zi = tp h (y), 

where O C Ai is a suitable open subset containing A4 ff . 

Definition 1.4 (Geometric RATTLE). Let 93^ be a method approximating exp(/iA"#) 
The RATTLE mapping R h : M H B z H> Z\ £ M H is defined implicitly by 

Z!€/C 5i nA^, 5i = 5 fc (z ). 

These are abstract definitions of SHAKE and RATTLE. In order to practically be able 
implement them, an implicit definition of M in terms of constraint functions, and a 
parameterisation of 1C ZQ , is needed. This issue is discussed in § 4, and is related to 
Assumption 2 introduced below. 

In the holonomic case it is already known that SHAKE and RATTLE essentially yield 
the same method, since the projection step at the end of RATTLE is "neutralised" by the 
projection step in SHAKE. This observation is made geometrically precise in § 4, where 
we show that SHAKE and RATTLE are two different representations of the same fiber 
mapping. 

Well-posedness of SHAKE and RATTLE The algebraic equations defining SHAKE and 
RATTLE can be thought of as discretisations of the original equation (1). However, 
contrary to the continuous case, the discretised equations are not intrinsic to A4. Thus, 
well-posedness of SHAKE and RATTLE depends on how A4 is embedded in V. 

Let TA4 denote the orthogonal complement of TA4 with respect to the symplectic 
form oj, i.e., u E T X A4 if and only if u(u, v) = for all v £ T X A4. 

Definition 1.5. A submanifold M of V is called coisotropic if TA4 1 - C TAf . 

As is explained carefully in §3, the natural assumption in order for SHAKE and RATTLE 
to be well-posed is the following, which is a completely extrinsic condition, i.e., it only 
has to do with how Ai is embedded in V . 



4 



Assumption 2 (Coisotropy). Ai is a coisotropic submanifold of V . 

Example 1.6. For the setting in Example 1.1, let gi, ■ ■ ■ ,g m be the components of the 
vector valued constraint function g. Then Assumption 2 means that 

{9i,9j}(z) = 0, Vz€«M. 

An equivalent interpretation of Assumption 2 is that none of the Lagrange multipliers 
in equation (2) are resolved by differentiating the constraint condition once. From a 
DAE point of view, the nondegeneracy and coisotropy assumptions together asserts that 
equation (2) has index 3. An important particular case is obtained when g does not 
depend on p. In that case, Assumption 2 always holds (see § 5.1). As shown in § 3.2, 
Assumption 2 also implies that fC z is parameterised by 

m 

(Ai,...,A 

m) 1 ^ 6Xp ( ^ ^ \, l Xg i J (z). 
i=l 

In turn, this means that the geometric SHAKE method is given by 

m 

S h = ip h o exp ( ^2 ^i x ai) 
i=i 

where Ai, . . . , \ m are determined implicitly by the conditions gi° Sh = 0. Likewise, the 
geometric RATTLE method is given by 

m 

R h = exp ( ^2 a i x ai) ° S h 
i=i 

where <ti, . . . , o m are determined implicitly by {Xjj ■ Vgi) o R h = 0. 

It is easy to find instances where the coisotropic and/or the nondegeneracy assump- 
tions do not hold, and where the SHAKE and RATTLE methods are not well-defined. 
For example, if we take as constraint H = const, then the nondegeneracy assumption 
does not hold, and it is easy to see that SHAKE and RATTLE are not well-defined. This 
is expected, since the result by Ge and Marsden [3] asserts that it is not (in general) 
possible to construct symplectic and energy preserving methods. In § 5.2 we give further 
examples of failing assumptions. Lastly, in § 5.4 we also give a numerical example of a 
Hamiltonian problem with mixed position and momentum constraints, where we use the 
geometric SHAKE and RATTLE methods. 

Main results The main results in the paper can be summarised as follows. 

1. Under Assumption 1, the set Ai H is a symplectic submanifold with symplectic form 
w = i\. H v, and equation (1) is well-posed for initial data in Ai H . (Theorem 2.7) 
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2. Under Assumption 1 and Assumption 2, there exists an open set O C Ai, contain- 
ing Ai H , such that the SHAKE map S^: O — > O is well defined and presymplectic, 
i.e., Sty = v. Further, it is convergent of order at least 1. (Theorem 4.1 and 
Proposition 4.4) 

3. Under Assumption 1 and Assumption 2, the RATTLE map Rh : Ai H — > Ai H is 
well defined and symplectic, i.e., R* h vo = w. Further, it is convergent of order at 
least 1. (Theorem 4.1 and Proposition 4.4) 

2 Hamiltonian Systems on Presymplectic Manifolds 

In this section we investigate the geometric structures of equation (1) from the intrinsic 
viewpoint, i.e., without "looking outside" of Ai. 

In general, a presymplectic manifold is a pair (AA,i>), where Ai is a smooth mani- 
fold, and v is a closed 2-form on Ai called a presymplectic form. The difference from a 
symplectic form is that v need not be non-degenerate. Thus, a symplectic manifold is a 
special case of a presymplectic manifold. We review some geometric concepts of presym- 
plectic manifolds that are essential in the remainder. For a more thorough treatment, 
we refer to the book by Libermann and Marie [12]. 

Given a function Hj^ £ C°°(Ai), equation (lb) constitutes a Hamiltonian system 
on (Ai,i>). Since v might be degenerate, this equation is not, in general, an ordinary 
differential equation, but instead a DAE. We show in § 2.3 that under Assumption 1 it is 
an index 1 DAE on Ai. (In § 3 we take the complementary extrinsic viewpoint, and we 
show that under Assumption 1 and Assumption 2 equation (1) can be interpreted as an 
index 3 problem on "P.) 

2.1 Foliation 

Throughout the paper we make the following "blanket assumption" : 



The dimension of the kernel distribution kerf is constant. 



One important consequence is that the distribution ker v (now assumed to be regular) 
is integrable (cf. [6, Th. 25.2]). That is, at each point x 6 Ai there is a submanifold 
KL X C Ai passing through x whose tangent spaces coincides with the distribution. The 
submanifolds fC x are called leaves, and the collection K, of all leaves is called a foliation. 

Remark 2.1. The foliation defines an equivalence class by y € [x] if y £ K x . We denote 
the set of all such equivalence classes by M. The projection is given by 

tt: M 3 x [x] e M. 

The set Ai may or may not be a smooth manifold. When it is, the presymplectic form v 
descends to a symplectic form v on Ai, and ir is a symplectic submersion. The projection 
map 7r being a submersion means that we have a fibration of Ai. Locally, every foliation 
is a fibration, but not necessarily globally. 
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A presymplectic mapping is a mapping ip: M —> M that preserves the presymplectic 
form v, i.e., for which 

v(u, v) = ip*u(u, v) := v(Tip ■ u, Tip ■ v) V u, v G T X M.. 

There is a certain class of mappings that are trivial in the sense that they reduce to 
the identity mapping in the quotient manifold Ai. 

Definition 2.2. A smooth mapping ip: Ai — > Ai is called trivially presymplectic if it 
preserves each leaf, i.e., if 

ip(x) G JC X for all x G M. 

The following result is clear. 
Proposition 2.3. If ip is trivially presymplectic, then it is presymplectic. 

2.2 Hidden Constraints 

The fact that v does not have full rank reflects that, in general, the possible solutions 
of (lb) do not fill the whole manifold Ai. Indeed, if a curve j(t) is a solution of (lb), 
then d-ff (7(t)) must be in the set v(Tj( t \Ai). As already seen in § 1, the set of points at 
which this is fulfilled defines the hidden constraint set Ai H C Ai, given by (3). 

Remark 2.4. In general, this set is defined as the locus of m functions, where m is the 
dimension of K, x . However, if the differential of those functions are not independent at 
the locus points, the set Ai H need not be a submanifold, and if it is a submanifold, it need 
not be of codimension m. For instance, we may have Ai H = M if the Hamiltonian H 
is constant along each leaf fC x . This is in particular the case if v is non-degenerate. 

Remark 2.5. The subset Ai H is, strictly speaking, not a set of hidden constraints, but 
rather implicit constraints as a consequence of (lb). 

2.3 Nondegeneracy Assumption 

In this section we show that the nondegeneracy assumption, Assumption 1, ensures 
that: (i) A4 H is a submanifold of A4; (ii) w = i* m h u is a symplectic form; and (iii) the 
initial value problem (lb) is a Hamiltonian problem on (Ai H ,w). As a consequence, 
problem (lb) have unique solutions for initial data in At . 

From a DAE point of view, Assumption 1 ensures that the DAE (lb) on Ai has index 1. 
As it turns out (see § 4 below), the nondegeneracy assumption, together with Assump- 
tion 2, also asserts that the geometrically defined SHAKE and RATTLE methods are well 
defined. 

We start with the observation that Ai H is the set of critical points of Hjc x - 
Proposition 2.6. For x G M, let Hjc x = i* Kx Hm- Then 

M H = {yeM:dH Kx (y) = 0} 
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Proof. If y G n /C x then dH Kx (y) = since dif(y) G i/(Tj,.M). Thus, the set 
Ai H n /Ca; consists of critical points of the function Hjc x . □ 

Theorem 2.7. Under Assumption 1, the following holds. 

1. The set A4 H is a submanifold of M. 

2. At a point x G At we /iai>e 

TjjAI = J x M H 8 ker v. 

In particular, the presymplectic form v restricted to A4 H is a symplectic form. 
Thus, Ai H is a symplectic manifold. 

3. Equation (lb) has unique solutions for initial data in A4 H . These solutions are 
given by the solutions of the Hamiltonian problem on M. H obtained by restricting H 
to M H . 

Proof. Each statement is proved respectively as follows. 

1. Let Xi, . . . ,X m be linearly independent vector fields on M that span the distri- 
bution kerzA Define the functions Pi(x) := (dH(x),Xi(x)). Then, using Proposi- 
tion 2.6, 

M H = { x G M : pi(x) = 0, i = l,...,m}. 

M. H is a submanifold if dpi(x), . . . , dp^(x) are linearly independent for every x G 
M, . An equivalent conditions is that the matrix 

niij := (dpi(x),Xj(x)) , i, j = 1, ...,m 

be invertible for every x G Ai H . Using that Pi(x) = 0, we get in local coordinates 
Xj that 

where Xi = Xf-J^. Since X\{x), . . . ,X^{x) is a linearly independent basis of 
(ker v) x = T x tC x , Assumption 1 means exactly that this matrix is invertible, which 
thus proves the first assertion. 

2. For the second assertion, using that the codimension of M is m, it suffices to 
prove that (kerz^ n T X M H = for every x G M H . Let u G T X M H . Then 
(dpi(x),u) = 0. Next, assume that U G (kerv) x . Then U can be expanded as 
U = u % Xi{x). We now get 

k 

= ^^"U 4 (dpj,Xi(x)) = miju\ 
i=i 

Under Assumption 1 we know that mij is invertible, which implies that U = 0. 
Thus, (ker v) x D T X A4 H = 0, which proves that v restricted to At is nondegener- 
ate. 
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3. For the final assertion, it is enough to show that j(t) is a solution to equation (lb) 
if and only if it is a solution to the Hamiltonian problem 

u(j(t),U) = {dH(i(t)),U) , VUeTM H 

on the symplectic manifold M H (for which existence and uniqueness follows from 
standard ODE theory). As we have seen, under Assumption 1 every X G T^^Ai 
can be written X = U + W with U G TM H and W G kerz/. Now, 



v(i(t),X) = i/( 7 (i), U) and (dH( 7 (t)),X) = (dff ( 7 (t)), tf> 
where the first and second equality follows respectively since W G ker v and 

(dH(rf(t)),w) = (dHfrWiYtV'xMt)) ) = ZX^M*)) = °- 



This ends the proof. □ 

Remark 2.8. If the prescribed initial condition does not lie in the set Ai H , there 
cannot be any solution curve passing through this point. On the other hand, if Ai H 
is a submanifold, and if it intersects the leaves K, cleanly, i.e., if the dimension of the 
intersection is constant, and if that dimension is larger than zero, then the equation may 
have infinitely many solutions. This is what happens if H is constant on the leaves of 
the foliation K. 

The following result will be useful in § 4, when we analyse SHAKE and RATTLE. 

Corollary 2.9. Under Assumption 1, there exists an open set O C A4 containing Ai H 
such that the equation y G K x n Ai H has a unique solution for every x G O. The 
corresponding trivially presyraplectic projection map IT: O — > Ai H , defined by n(rr) = y, 
is a submersion. 

Proof. This follows from Theorem 2.7 item 2, namely that for x G Ai H , T x Ai = 
T X M H ekerz;. □ 



3 Coisotropic Constraints 

In this section we study the geometry of problem (1) from the extrinsic viewpoint. That 
is, we study properties of M. as a submanifold of the symplectic manifold (V,uj). Notice 
that v := i* m oj is a presymplectic form on A4, since du = di* M to = i* M dui = 0. Thus, any 
submanifold of a symplectic manifold is automatically a presymplectic manifold. 
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3.1 Lagrange Multipliers 



Typically, a constraint manifold is defined in terms of a number of constraint functions. 
To this extent, let V be a vector space of dimension m, and denote by V* its dual. Let 
J : V — > V* be a smooth function such that the constraint submanifold Ai is given by 

M = r 1 (0) = {zeV:S(z) = 0}. (4) 

If is a regular value for J, i.e., if TJ(z) has full rank for all z 6 V such that J(z) = 0, 
then A4 is indeed a regular submanifold of V . The dimension m of V is the number of 
constraints, i.e., the codimension of M. 

The problem (1) may now be reformulated as finding a smooth curve 

tt-> (*(*), A(t)) eV xV 

such that 

u,(i) = dtf(z) + d(J,A>(z), 

(5a) 

= J(z). 1 ; 

Here, the notation (J, A) means the smooth function z i— > (J(z),A), depending on the 
parameter A. The equation can equivalently be written as 

z = X H (z) + X {lA) (z), 

(5b) 

= J(*). 

We sometimes single out a basis {ej}j = i r .. im of V and define the functions gi by 

gi(z) := (J(»,ei) . (6) 

Notice that in the case T 3 = R , equation (5) coincides with equation (2) in Example 1.1 
above, with A = (Ai, . . . , A m ) being the coordinate vector of A, i.e., A = Y%Li ^« e «- 

The system (5) is again a DAE. Under Assumption 1, it follows from Theorem 2.7 
above that this DAE has unique solutions for initial data in Ai H . From a DAE point of 
view, Assumption 1 asserts that system (5) has index 3. 



3.2 Coisotropy Assumption 

Due to the solvability result imposed by Assumption 1, the Lagrange multipliers may 
be resolved as functions of z, which turns equation (5) into 

m 

z = X H (z) + ^ \i(z)X gi (z) =: X(z). (7) 
i=i 

Notice that X{z) is only defined for z £ M H and also that X(M H ) G JM H , so X(z) 
defines an ODE on the hidden constraint manifold At . From Theorem 2.7 it follows 
that its flow is symplectic. However, the individual vector fields \i{z)X gi {z) are not 
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Hamiltonian vector fields on V (assuming that \{z) is denned also outside of Ai ). In 
this section we present an assumption on the embedding 1^4 : Ai — > V which ensures that 
vector fields of the form f(z)X gi (z) are trivially presymplectic vector fields on Ai. As 
we will see in § 4, this is essential in order to ensure presymplecticity and symplecticity 
of SHAKE and RATTLE. 

Recall from Definition 1.5 that Ai is a coisotropic submanifold of V if TA4 -1- C TAi. 
Also recall Assumption 2 above (the coisotropy assumption), which states that At is a 
coisotropic submanifold. We continue with some consequences of Assumption 2, which 
are later used in the geometric analysis of SHAKE and RATTLE. 

Remark 3.1. It is straightforward to verify that Ai being coisotropic is equivalent to 
TAt being isotropic, i.e., such that uj restricted to T At 1 - is zero. 

Remark 3.2. In a sense, a coisotropic submanifold is such that the symplectic form 
becomes as degenerate as possible (given a fixed number of constraints) when restricted 
on the submanifold. More precisely, a coisotropic submanifold is such that the dimension 
of the distribution kerz^, i.e., dimension of the leaves of JC, is equal to the number of 
constraints m. 

Remark 3.3. Prom a theoretical point of view, Assumption 2 is not a restriction on the 
presymplectic manifold Ai , since every presymplectic submanifold may be coisotropically 
embedded in a symplectic manifold [5]. 

Remark 3.4. In practice as shown in Proposition 3.5, a sufficient condition for the 
manifold Ai defined by the equations gi(z) = for 1 < i < m to be coisotropic is simply 
that 

{9i,9j}{z) =0, ij = l,...,m, \/ z e Ai. 

In particular, if the manifold Ai is defined by one constraint, i.e. if m = 1, then it is 
automatically a coisotropic submanifold. 

The following result gives alternative characterisations of coisotropic submanifolds. 

Proposition 3.5. Suppose that Ai is a submanifold ofV . Then the following conditions 
are equivalent. 

1. Ai is a coisotropic submanifold, i.e., TAi^ C TAi. 

2. keriy = JM L 

Further, if Ai is defined implicitly by (4), then the conditions are also equivalent to 

3. For any a, (3 £ V, the functions (J, a) and (J,/3) are in involution on Ai, i.e., 

{(J,a),(J,/3)}(z)=0 for z G At. 

4- For any a £ V , the Hamiltonian vector field Xn^ a \ is tangent to Ai. 
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In order to prove this, let us start with a lemma concerning the span of the Hamiltonian 
vector fields ^(j, a )- 

Lemma 3.6. Define the distribution 

O = span{X (J)a> (.M) : a G V}. 

Then JAi 1 = O. 

Proof. We show that O 1 - = TAi, which is equivalent to the claim. 

X G JM <=^ (d(J,a) ,X) = VaGF ^ w(I (Ji(l) ,I) = Va G V 

□ 

Proof of Proposition 3.5. We do it step by step. 
102 In general, 

kerz/ = JM n T«M J ", 

so kerf = T.M T.M 1 - C JM, and that is the definition of coisotropicity 

of M 

l-H-3 First, for x G M 

{gi,9j}(x) = <^ = 0, 

so the functions gi are in involution on M if and only if O (defined in Lemma 3.6) 
is isotropic, which is equivalent to M being coisotropic. 

3o4 Finally, it suffices to observe that for a point x G M, 

x {la) (x)eJ x M ^ (d(J,p),x {la) )(x) = o \f(3ev 

{(J,a),(J,/3)}(x) = V/3 G V. 

□ 

The following results follows directly from Lemma 3.6 and Proposition 3.5. 
Corollary 3.7. Let x G M. Then, under Assumption 2, the map 

V B exp(X( J)a ))(x) G K x 
is a local diffeomorphism. The fibre IC X is thus locally parametrized by V . 
Corollary 3.8. Let f G C 00 ^) and a £ V. Under Assumption 2 the vector field 

X(z) := f(z)X {la) (z) 
is tangent to M, and presymplectic when restricted to Ai. 
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3.3 Relation between X H and M H 

As a subset of V, the hidden constraint set A4 H is given by the points on A4 where the 
Hamiltonian vector field Xh is tangential to Ai. Let II denote the projection onto Ai H 
defined in Corollary 2.9. 

Proposition 3.9. Under Assumption 2, the hidden constraint set Ai H is 

M H = {zeM: X H (z) G J Z M }. 
Moreover, the differential equation (7) on Ai H can be written 

z = J z U-X H (z). 

Proof. 1. If z G Ai H , then by definition there exists Y G T Z AA such that 

(dH,X) =u(Y,X) VAgT z AL 
Since u(Xjj) = dH, that is equivalent to 

Xh - 76 TA^" 1 . 

Noticing that Assumption 2 means that TA4 1 " C TA^, and using Y G TA4 yields 
X H G T,A^. 

2. The differential equation on A^^ is such that 

A) = (dff, A) VA G J Z M 

so we obtain 

z-Ye TzM 1 , 

and z = T Z H ■ Xh(z). 

□ 

Remark 3.10. There are now several ways to compute Ai H . First, without any as- 
sumption, one can use the definition (3), and its immediate consequence Proposition 2.6. 
Under Assumption 2, one can also use Proposition 3.9. If the constraint manifold A4 is 
defined as in (4), a further useful description of M. H is 

M H = {z G M : {g h H}(z) = i = l,...,m}. 
This follows from the observation that 

X H G T,A^ <=> (d gi ,X H ) =0, i = l,...,m (8) 

and (dgi,X H ) = {gi,H}. 

Based on Theorem 2.7, Proposition 3.5 and Lemma 3.6, we can say much more on the 
behaviour of Xh in a neighbourhood of Ai H . Indeed, we have the following result, which 
is a key ingredient in the well-posedness of SHAKE and RATTLE, as will be explained in §4. 
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Lemma 3.11. Let y £ M and define the function 

F: K y — > V* 

x i — > dJ(x) • Xh(x). 

Then, under Assumption 1 and Assumption 2, the differential of F at y, i.e., the linear 
mapping 

dF(y) : J y K y — > V* , 

is invertible. 

Proof. In terms of the previously introduced basis {a}i = i m , the function F is given 
by 

m m 

F(x) = (d gi (x),X H (x)) e t = -J2 (dH(x),X g .(x)) e, 
i=i i=i 

Under Assumption 2, it follows from Proposition 3.5 and Lemma 3.6 that X gi (y) , . . . , X 9m (y) 
is a basis for T y K y . Relative to this basis, and the basis {ei}i=i s „. jm of V, the Jacobian 
matrix of dF(y) is given by 

ma ■= (d(dH(y),X gi (y)),X g .(y)) = {{H, 9i }, g 3 }(y) 

Define pi{x) := (dH(x), X gi (x)). Then rriij = (dpi{y),X g . (y)). Since y G Ai H we have 
that Pi(y) = 0. Now, under Assumption 1 and the exact same argument as in the proof 
of Theorem 2.7, it follows that rriij is invertible. This concludes the proof. □ 

4 Geometry of SHAKE and RATTLE 

Geometrically, the basic principle of SHAKE, defined in Definition 1.3, with ip^ as an 
underlying method, can be described as follows. For some initial data z € At, slide 
along the fibre with z + £ KL Z such that (ph(z + ) "lands" again on the submanifold A4. 
The RATTLE method, defined in Definition 1.4, is then a post-processed version of SHAKE, 
which is described geometrically as follows. For some initial data zo £ A4 H , take one 

step with SHAKE landing on z7 £ M.. then slide along the fiber JC - to end up on z\ £ 

z i 

A4 H n /C -. This process is visualised in Figure 1. 

z i 

If we assume that the SHAKE map Sh '■ O — > O is well defined for some open subset 
O C A4 containing JV[ H , then we may define a sliding map T iph = (pT oSh- Consequently, 
it follows from Definition 1.3 that the sliding map O 3 z \— > T iph (z) £ O is defined 
implicitly by the equation 

Mz + )^M, z+£/c 2 no, r (ph (z) = z+. (9) 

Notice that: (i) T !fih is fiber preserving, i.e., presymplectic, and (ii) T lfih is a projection, 
i.e., T Vh oT lfih = T Vh . Since Sh = ^Ph ^f h it follows that SHAKE, if it is well-defined, is a 
fiber mapping, i.e., it maps fibres to fibres. It is, in fact, a little bit more than that, since 
it maps the whole fibre tC z n O to the same point Sh{z). Hence, when using SHAKE it is 
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Figure 1: An illustration of the SHAKE and RATTLE method. The first point z must lie in 
the vicinity oi M. H . The subsequent points z± , z^ , ■ ■ ■ produced by SHAKE will 
all lie in the modified hidden constraint manifold Ai H . The points zi, z%, . . . 
produced by RATTLE will lie on At . Notice how the points z7, z^ and z^ 
always lie on the same leaf. Moreover, the precise location of a point on the 
fibre (for instance z on the figure) is irrelevant for both SHAKE and RATTLE. 
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not important where on the initial fibre K z one starts (as long as it is close enough to 
A4 H so that SHAKE is well-defined). Furthermore, regardless of where on the fiber one 
starts, after one step SHAKE remains on the modified hidden constraint set, given by 



Since T Vh is a projection, A4 is strictly smaller than O. If SHAKE is well-defined, M 
is in fact a symplectic submanifold of A4 (see Proposition 4.7 below). 

Let II be the projection on M H given in Corollary 2.9. Assume that II is well-defined 
on O. Then RATTLE is given by Rh = II o Sh- Notice that SHAKE and RATTLE defined 
exactly the same fiber mapping. In particular, 



Thus, RATTLE is only a cosmetic improvement of SHAKE, and has no influence on the 
numerical scheme except at the last step. 

We now give explicit conditions under which SHAKE and RATTLE are well defined and 
can be computed. More precisely: 

1. When is a method iph such that the corresponding SHAKE and RATTLE methods 
are well defined? 

2. How can we parameterise Y Vh (so that Sh is computable)? 

3. Will SHAKE and RATTLE converge to the solution of equation (1) as h — > ? 

4. Are Sh and Rh presymplectic as mappings O — > O ? 

5. Is Rh symplectic as a mapping M. H — >■ M. H , and Sh symplectic as a mapping 
M H ^M H ? 

6. Can Sh or Rh be reversible? 

These questions are addressed in the remainder of this section. 

4.1 Well-Posedness 

In order for SHAKE and RATTLE to be well defined, we need the "sliding process" to 
have a locally unique solution. Whether so or not depends on the map iph- 

Theorem 4.1. Suppose that Assumption 1 and Assumption 2 hold. Consider a (smooth) 
method iph, consistent with z = Xjj(z). Then for h small enough and for z G M. in a 
neighbourhood of M. H the equation 



(f h oT Vh {0). 



Uo(S h ) k =Uo(R h ) k . 



(10) 



j(^(z+))=o, z + eic z 



(ii) 



has a unique solution. 



16 



Proof. Let x G M, . We define the function Fh'- K, x —¥ V* by 







ah 

We see that Fh depends smoothly on h. Notice that for h ^ 0, 

J(^(*+))-% (* + )) J(^ + )) 



h h 



because tpo(z + ) = z + and z + G M. J(z + ) = 0. 

Consider now the case h = 0. The method (fh is consistent, so (d9?h/d/i)|/,, = o(.z + ) = 
Xh(z + ), which shows that 

F (z + ) = dS{z + )-X H {z + ). (12) 

We want to find a neighbourhood Oh C IC X oi x such that the equation Fh{z + ) = 
with z + G has a unique solution for small enough h. This will prove the claim. 
The strategy is to show that dFh(x): J X 1C X — > V* is non-singular. We start with the 
case h = 0. 

Using (12) and appealing to Lemma 3.11, dFo(x) is invertible under Assumption 1 and 
Assumption 2. Thus, by the inverse function theorem, we can find an open neighbour- 
hood Oq C K, x such that Fq: Oq — > Fq(Oo) is a diffeomorphism. Also, since F^ depends 
smoothly on h, it follows that dFh(x) is invertible for small enough h. Thus, for small 
enough h, we can find an open neighbourhood Oh C K, x such that Fh ■ Oh — > Fh(Oh) is 
a diffeomorphism. 

It remains now to show that G Fh(Oh), so that the equation Fh(z + ) = with 
z + G Oh has a unique solution. Now, if x G Mr it follows from Proposition 3.9 that 
Xh(x) is tangential to At, which means that dJ(x) • Xh(x) = 0, so we get Fq(x) = 0. 
Thus, G Fq(Oo), and it follows by smoothness that G Fh(Oh) for small enough h. □ 

Remark 4.2. The coisotropy assumption is essential for the result of Theorem 4.1, 
because it uses Lemma 3.11 which depends on that assumption in an essential manner. 
On the other hand, if (11) has a unique solution in a neighbourhood of Ai H , then 
Assumption 2 must hold. 

Indeed, the theorem implies that K, x is locally diffeomorphic to V*. This implies that 
their dimension is the same, so ker v must have the same dimension as V* . 

As a result, 

dim ker v = d'miV* = codimA'f = dimT.A/f"'". 

Now, in general ker v = TM D TAi 1 ', so we get TM 1 - C TA4, which implies that M. is 
a coisotropic submanifold. 

The result of Theorem 4.1 allows to define the sliding map T lfih by (9). The corre- 
sponding SHAKE and RATTLE maps, Sh and Rh, are thus defined when h is sufficiently 
small. 
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4.2 Fibre Parametrisation 

When the manifold Ai is defined implicitly as the locus of functions gi, one can express 
the effect of T Vh using the flows of the Hamiltonian vector fields X gi . 

Proposition 4.3. Suppose that the p^-SHAKE method is well defined, so that T Vh is 
well defined. Under Assumption 2 there exists functions Ai, . . . , A m G C°°(.M) such that 
the sliding map T Vh is given by 

m 

T Vh (x) = exp (h -^0)^) (x) 
i=i 

Proof. Follows directly from Corollary 3.7. □ 



4.3 Convergence 



We may now show the convergence of SHAKE and RATTLE. The proof is essentially 
a standard convergence argument for RATTLE, which is a numerical method on the 
manifold A4 H . The convergence of SHAKE is then obtained using (10). 

Proposition 4.4. Suppose that Assumption 1 and Assumption 2 hold. Let ipi t be a 
method consistent with z = Xh{z). Then the phSHAKE and ip^-RATTLE methods are 
convergent of order at least 1 . 

Proof. We first turn to RATTLE. The continuous system is an ordinary differential equa- 
tion on A4 H with vector field Til o Xh (see Proposition 3.9). Since it is an ordinary 
differential equation, we only need to show that RATTLE (defined in Definition 1.4) is 
consistent and standard arguments may then be used to show convergence of order 1 
(see [7, § II.3]). Using that Rh = II o p h o T iph and differentiating at h = we get 



dR h 



dh 



h=0 



Til 



dipt, 



dh 



dr 



•rh 



h=0 



dh 



h=0 



which follows since po = id and 

dR h 



dh 



id. The second term is in the kernel of Til, so 

diph 



h=0 



TLT 



dh 



h=0 



The assumption that ph is consistent with z = Xh(z) means that ^ L \ h _ = Xh- 
Consistency of Rh then follows. 

Now, using (10) and that T iph (z) = z + 0(h) when h — > 0, we also get that SHAKE 
converges of order at least 1 . □ 

Remark 4.5. Interestingly, the local error of SHAKE is only 0(h), so standard arguments 
do not apply to show convergence. However, due to the fact that SHAKE is the same fibre 
mapping as RATTLE, this error does not accumulate, and the global error is still 0(h). 
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4.4 Symplecticity 



We examine in which sense SHAKE and RATTLE may be regarded as symplectic method. 
The essential result is that both SHAKE and RATTLE are presymplectic, i.e., they preserve 
the presymplectic structure of M.. 

Proposition 4.6. Let iph be a symplectic method. Then the corresponding SHAKE map 
Sh and RATTLE map Rh, regarded as mappings M. — > M, are presymplectic. 

Proof. Since T lfih is trivially presymplectic, it is in particular presymplectic, so 

u(u, v) = v(Tr,p h ■ u, TT Vh -v), V u, v £ JM. 

Further, since ifh is symplectic it follows that 

v(u, v) = u(Jip h • (Tr^ • u), l(p h ■ (JV Vh ■ v)) 

= v{J{ip h o r^j • u, T(tp h o r^j • v) 

= v (TSh ■ u, TSh • v) , V u, v £ TM. . 

Thus, Sh is presymplectic. 

Moreover, since Rh = II o 5^ and IT is trivially presymplectic, Rh is also presymplectic. 

□ 

Proposition 4.7. Let ifh be a symplectic method. Under Assumption 1 and Assump- 
tion 2, the set Ai H is a symplectic submanifold of M., with symplectic form vb = l *j^ H 1 '- 

Proof. Under Assumption 1 it follows from Theorem 2.7 that the set (Ai H ,w) is a 
symplectic submanifold of M. We first show that M is diffeomorphic to M. , and 
thus also a submanifold of A4, by constructing a diffeomorphism M H — > Ai H . Under 
Assumption 1 and Assumption 2, the map T iph : Ai H — > A4 is well defined. Let z £ A4 H . 
By construction it holds that ker(T z T Vh ) = T Z JC Z . Thus, J z T iph : T Z M H —> Tp^ ( Z )M 
is injective, so T Vh : Ai H — y T iph (A4 H ) is a diffeomorphism. Next, since ifh - V — >■ V 
is a diffeomorphism, it also holds that ifh - T, Ph (M H ) — >■ iph(T [ph (M. H )) = M H is a 
diffeomorphism. Thus, Sh = <fh ° ^<p h is a diffeomorphism as a map Ai H — > M. , so 
M H and M H are diffeomorphic. 

Next, we show that the form w = l *j^ H v is non-degenerate. Let y = Sh{z) and 

u, v £ T y A4 H . From Proposition 4.6 it follows that Sh is presymplectic, so 

V(U, V) = ^(TyS^ 1 ■ U, TyS^ 1 ■ v) = W^yS^ 1 ■ U, TyS^ 1 • v) . 

Since TyS^ 1 : Ai H — > M. H is invertible, and since (Ai H , w) is a symplectic submanifold, 
it follows that v(u, v) = for all v £ T y Ai H only if u = 0, which shows that (A4 H ,z&) 
is a symplectic manifold. □ 

The following result follows directly from Theorem 2.7, Proposition 4.6 and Proposi- 
tion 4.7. 
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Corollary 4.8. Let ip^ be a symplectic method and assume that Assumption 1 and 
Assumption 2 hold. Then: 

• The SHAKE map Sh, regarded a mapping A4 H — > M , is symplectic. 

• The RATTLE map Rh, regarded as a mapping At —> M , is symplectic. 

4.5 Time Reversibility 

Just as in the holonomic case, if the underlying method tp h is symmetric, i.e., ip7 = ip~h, 
then RATTLE is also symmetric and thus of second order. Note that SHAKE can never 
be symmetric, because although Sh preserves Ai H , the reverse method S-h does not. 

Proposition 4.9. If the underlying method <ph is symmetric, then so is RATTLE, con- 
sidered as a method from A4 H to Ai H . 

Proof. This follows from the symmetry property of the map Tm h . In general, if z^ = 
T v {zq), zT = ip(zq), and Z\ = H-(zT) (see Figure 1), then 

z7 = r^-i^i). (13) 

This follows from Theorem 4.1, because z^ is a solution of the equation j((/3 -1 (zi)) = 0, 
so it must be equal to T ip -i(zi). 

Suppose that we start from a point zq G A4 H . The image by the RATTLE map is 
z\ = Rh(zo). Now, since zq = H(zq ) (as we assumed that zq G M h ), and (13), we 
obtain zq = H o tp^ 1 o r^-i. 

If we assume now that (fh is symmetric, i.e., ipj^ 1 = tp~h, we obtain zq = R^h(zi), so 
RATTLE is symmetric. □ 

5 Examples 

In this section we give examples of constrained problems that can be solved with the 
geometric SHAKE and RATTLE methods. In §5.2 we also give non-trivial examples where 
the non-degeneracy assumption fail to hold. 

5.1 Holonomic Case 

In this section we study the classical, so-called holonomic case, where the constraints 
depend only on the position q, and not on the momentum p. 

Consider the classical case of constrained mechanical systems, where the symplectic 
manifold V is simply R 2n , with coordinates (q l ,Pi), equipped with the canonical sym- 
plectic form u) = 'Yui&q 1 A dpi. The constraints are given by g(q) = for a function g 
defined from V to V := R m . It should be emphasized that g does not depend on p, 
i.e., g p = 0. The leaf passing through the point z = (q,p) is then an affine subspace 
described parametrically by 
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Assumption 2 is automatically fulfilled, because ^ = implies 

{9i,9j} = 0. 

Moreover, the standard assumption placed on the Hamiltonian H is that the matrix 
g q (q)Hpp(q,p)g q (q) T be invertible, (14) 
(see [8, § VII. 1.13], [15]), which implies Assumption 1. 

Remark 5.1. The condition (14) is in fact more stringent than Assumption 1, since 
we only need it at critical points of H on a leaf, i.e., on the points lying on the hidden 
constraint manifold Ai H . It is clear that the condition (14) is too restrictive when the 
leaves are compact. We will examine such an example in § 5.4. 

5.1.1 Classical SHAKE and RATTLE 

The classical SHAKE method is only defined for separable Hamiltonian functions (see [8, 
§ VII. 1.23], [15], [11]). However, it is readily extended to general Hamiltonian as the 
mapping (q ,Po) ^ (QIiPi) defined by 

Po = Po - ^g q (qo) T ^ 

Pi/2 = Po - ^H q (q ,p 1/2 ) 

qi = qo + ^{H p {qo,Pi/ 2 ) + H p (q 1 ,p 1/2 )) 

Pi = Pl/2 ~ ^H q (q 1 ,p 1/2 ) 
= g(qi). 

The classical RATTLE method is obtained by adding a projection step onto the man- 
ifold M. H , i.e., the next step is instead (qi,p±) i-> {qi,Pi) defined, on top of SHAKE, 
as 

Pi =Pi +g q (qi) T v 
= 9q(qi)H p (qi,pi). 

We see that the mapping ((ftjjPo) l— ^ {QIiPi) * s the Stormer-Verlet method, so the 
classical SHAKE and RATTLE methods are obtained using the Stormer-Verlet method as 
the underlying unconstrained symplectic method. 

Similarly, in [8, § VII. 1.3], a first order method is described as 

Po ~ hg q (q ) T X 
Po ~ hH q (q ,Pi) 
q + hH p (q 0l p-) 

g(qi) 




qi = 

= 
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We see that the mapping (go>£>o~) l— ^ (qi,P\) is the symplectic Euler method, so this 
is the SHAKE method obtained using symplectic Euler as the underlying unconstrained 
symplectic method. When complemented with the projection step 

Pl = p - - hg q (qi) J 'n 
= 9q(qi)H p (qi,pi), 

that method becomes the corresponding RATTLE method. 

5.2 Examples of Failing Nondegeneracy Assumption 

Let us examine an example where Assumption 1 is not fulfilled. The manifold V is R 4 
with coordinates q,p,q,p, and the symplectic form is w = dq A dp + dq A dp. The 
constraint function J is 

■)(q,p,q,p) = q- 

The corresponding leaves are the lines parametrised by p, i.e., the leaves are the lines of 
equation (q,p,q) = (qo,Po, %)• We choose the Hamiltonian 

H = — +pq- 

The hidden constraint manifold is then A4 = {q = }. The restriction of H on a leaf 
is a linear function (because H is linear in p), so its critical points are degenerate, and 
Assumption 1 is not fulfilled. As a result, there are extra constraints which prevent the 
problem from being well-posed on Ai H . It is readily verified that the system has in fact 
index five instead of three, the tertiary and quaternary constraints being respectively 
p = and p = 0. 

If we had assumed instead that H = the hidden constraint set would be Ai H = M, 
because H is now constant on the leaves (because H does not depend on p). In that 
case, there are infinitely many solutions to the problem at hand. 

5.3 Index 1 Constraints 

Let V = R 2n + 2fc with coordinates q G R n , p G R n , a G R fc , /3 G R fc , symplectic form 
uj = dq A dp + da A d/3, Hamiltonian H(q,p, (3), and consider the holonomic constraint 
a = 0. Because the constraints are holonomic, the coisotropy assumption is satisfied 
(see § 5.1). The presymplectic manifold is Ai = R, 2n + fc with coordinates (q,p,/3) and 
presymplectic form v = dq A dp. The hidden constraint submanifold is A4 H = { (q,p, A) : 
Hp(q,p,/3) = 0}, giving the presymplectic system 

q = H p (q,p,P) 

p = -H g (q,p,f3) (15) 
= Hp(q,p,P) 

The nondegeneracy assumption is that Hpp is nonsingular on At . When this holds, 
the secondary constraint can be solved for (3 = B(q,p). This is the case of index one 
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constraints considered in [14]. In this paper, a symplectic integrator is introduced and 
studied for the index one constrained Hamiltonian system (15), consisting of a direct 
application of a symplectic Runge-Kutta method to the full system (15). In the case 
when cph is the midpoint rule, we now show that this method is equivalent to the appli- 
cation of SHAKE or RATTLE applied on V. Indeed, Hamilton's equations on V are given 

by 

q = H p (q,p,p) 
P = -H q (q,p,P) 
a = Hp(q,p,P) 

p = o 

for which the midpoint discretization, defining iph, is 



gl 


- qo 
h 


= H P (- 


Pi 


-Po 

h 


= ~H q 






= Hp{ 




h 


Pi 


-Po 


= 



go + qi po+pi Po + Pi 
2 ' 2 ' 2 
go + gl Po + Pi Po + Pi 



(16) 



while the constraint flow is given by exp(A T X a ) : (q,p, a, P) {q,p, a, P — A). Initially, 
the primary constraints are satisfied, i.e. ao = 0, and the constraint flow determines Pi 
so that ax = 0, i.e., so that 

= H ( qo + qi Po+Pi Al+AA 
2 ' 2 ' 2 / 

with solution = i3( g °^ gl , Po t pi ) . This is exactly the method of [14] in the case of 

the midpoint rule. The final step of rattle chooses P so that = H^(qi,pi, P), but (as 
we have seen in the general case), this value is irrelevant for the next step. 

5.4 Harmonic Constraints 

Consider the symplectic manifold V defined by 

V := C 2 ~ T*R 2 , (17) 

with coordinates 

^o = go + ipo, «i = 9i + ipi, 
and the canonical symplectic form 

w(go,Po,gi,Pi) = dgo A dp + dgi A dpi- 
We will sometimes use the notation q = (qo, qi) and p = (po,pi). 
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The constraint space is V = Ft, and the single constraint function is 

s(q>p) : = Nil 2 + IIpII 2 - 1 - 

Note that the corresponding manifold Ai is S 3 , i.e., the unit sphere in R 4 . It is 
automatically a coisotropic submanifold since it has codimension one (see Remark 3.4). 
The foliation induced by g is the Hopf fibration of the 3-sphere in circles. The Hopf map 

(zo, zx) (->■ (IzqZi, \z \ 2 - \z x \ 2 ) (18) 

projects C 2 onto R 3 = C x R and maps 3-sphere onto 2-spheres. The quotient manifold 
Ai is thus diffeomorphic to a 2-sphere. 

The flow of the Hamiltonian vector field X g is 

ex.p(X g t)(z , zx) = (e^ce"^). (19) 
5.4.1 Hamiltonian without Potential Energy 

In this section, we examine the validity of Assumption 1 for the simple Hamiltonian 
function given by 



#(q,p) = 




The equations (5) are in this case 

q=(l + A)p 

p = -Aq (20) 
i II ii 2 , n i|2 

i = Nl + IIpII ■ 

Let us look at the restriction % of H on a leaf parametrized by 6 £ R/Z, passing by 
the point (q, p) for 6 = 0. We obtain after trigonometric simplification 

n(0) = ^(||P|| 2 (1 +cos(20)) + ||q|| 2 (l-cos(2^)) + 2q • psin(20)) . 

One checks that V,(0) = -ff(q, p) = ||p|| 2 /2. By differentiating, one obtains 

n'{6) = \ ((||q|| 2 - ||p|| 2 ) sin(20) + 2q • p cos(20)) . (21) 

We know that the point (q, p) lies in Ai H if 1~L'(0) = 0, which leads to 

M H = {(q ,po,q 1 , Pl ) eM :q-p = 0}. (22) 

This could also have been derived by differentiating the constraints in the differential 
algebraic equation (20), or by using one of the equivalent formulations of Remark 3.10. 

Such a point is a nondegenerate critical point for % if and only if ||q|| 2 — ||p|| 2 ^ 0. 
We also see that if (q, p) is a degenerate critical point, then q • p = and ||q|| = ||p||, 
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Figure 2: An illustration of the intersection of a leaf with M . The two opposite points 
on each leaf are part of the same connected component of Ai H . 

so H is constant on the whole leaf. A closer examination shows that there are two such 
leaves, which are mapped to the 2-sphere to the antipodal points (i, 0) and (— i, 0) by 
the Hopf map. We call those two points the exceptional points. 

Let us turn to the intersection of Ai H on each leaf. The form of (21) shows that % 
has four distinct critical points on each leaf. Moreover, two such points are such that 
||q|| > ||p|| where as the two others are such that the opposite inequality holds. The 
points are intertwined, see Figure 2. Interestingly, it is possible to join two opposite 
points on the fibre by a homotopy curve of the following form. If ||q|| < ||p||, then we 
choose the homotopy 

[-1,1] 3t^ (tq,p), 

whereas if ||p|| < ||q||, we choose the obvious other homotopy, where p and q are 
swapped. This shows that if one removes the leaves above the two exceptional points, 
the set Ai H is a manifold with two connected components. 

Note that ||q|| and ||p|| are constant on each trajectory, so the trajectories never meet 
the exceptional points. 

Let us gather our findings: 

Proposition 5.2. Consider the symplectic manifold V defined by (17), and the con- 
straint submanifold defined by 

M := { (zo,zi) G V : |z | 2 + N 2 = 1 and {2z z 1: |z | 2 - \z x \ 2 ) / (±i, 0) }. 

Then both Assumption 1 and Assumption 2 hold. The hidden constraint manifold Ai H 
is defined by (22) . It has two connected components. Each of those components crosses 
each leaf of K, twice, as depicted on Figure 2. Moreover, the vector field defined on Ai H 
is complete, i.e., all the trajectories exist for all time. 

5.4.2 Hamiltonian with Linear Potential Energy 

We now choose the Hamiltonian 

#=^-g-q (23) 

where g is a fixed given vector. 
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(a) Energy plot (b) Trajectories 

Figure 3: Some simulation results for the Hamiltonian (23) with the vector g := (0, — 1). 

Both the SHAKE and RATTLE methods use the midpoint rule as underlying un- 
constrained symplectic integrator. The time step is h = 0.1 in all the examples. 
We choose the initial conditions in Table 1. The plot (a) shows an energy plot 
for the SHAKE and RATTLE methods for the initial condition zq. The plot (b) 
shows the three trajectories for the initial conditions zq, z\ and Z2- We first 
use the Hopf map described in (18) and then use a stereographic projection. 

Now the structure of the hidden constraint manifold Ai H is more involved. In par- 
ticular, there are many singular points where the set Ai H ceases to be a manifold. For 
that reason it is somewhat difficult to find entire trajectories which stay on the patches 
of At which are manifolds. We nevertheless present in Figure 3 simulated trajectories 
for initial conditions such that those trajectories do not meet any singular point. 

6 Discussion 

The geometric version of the Dirac constraint algorithm [4] in the instance used here 
delivers a chain of submanifolds 

M H «-»■ M ^ V 

in which Ai H is symplectic. The geometric rattle algorithm delivers symplectic inte- 
grators on A4 H . However, V, the intervening presymplectic manifold Ai, its coisotropy, 
and knowledge of its fibers, are essential to the algorithm. If the fibers cannot be explic- 
itly parametrized, the algorithm is still formally defined, but more computation would 
be required in practice — for example, by integrating the fibers numerically to round- 
off error. This is an extreme version of a situation common in numerical analysis, in 
which allowing a wider class of methods (e.g. implicit Runge-Kutta methods, for which 
implicit equations have to be solved numerically) enables a wider class of properties. 
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zo 




Z3 


-0.7865261200000000 


-0.4963624948824013 


0.3477491188213400 


-0.4043988000000000 


-0.7319740436366664 


-0.8131619010029159 


-0.3880746864163783 


-0.4275775933953260 


-0.4368285559113795 


0.2173391755798215 


0.1225384882604160 


-0.0837800227934176 


-1.3703247300000001 


-1.3703247300000001 


-1.3703247300000001 



Table 1: Three initial conditions used in Figure 3. The fifth component is the Lagrange 
multiplier. 

If Ai is not coisotropic, then the coisotropic embedding theorem [5] says that there 
exists a symplectic manifold V' such that Ai is coisotropically embedded in V' (Re- 
mark 3.3). Thus, abstractly at least, one can extend the Hamiltonian on Ai arbitrarily 
to V' and apply the geometric rattle algorithm, for the rest of the required structure 
is instrinsic to Ai. In specific examples it may be possible to carry this out by finding 
a suitable symplectic vector space V' . The same remark holds if the given data is a 
Hamiltonian on a presymplectic manifold, see § 5.3. 

However there remain many constrained problems which do not fall into the classes 
considered here. The most fundamental one has data (V, Ai' , H) where Ai 1 is a symplec- 
tic submanifold of the symplectic manifold V . We do not know of symplectic integrators 
for this problem. They would provide symplectic integrators for a wide class of sym- 
plectic manifolds. A very general situation is that provided by the geometric version of 
the Dirac constraint algorithm [4], which, from presymplectic data (Ai, oj,H), produces 
a nested sequence of submanifolds 

Mr ^ Mr-i ^-t Mi := Ai 

defined by 

Ai k+i ■= {x £ A4k '■ dH(x) £ u(AAk) } 

and a (possibly nonunique) vector field X such that %m k {ix& — dH) = 0. One would like 
to integrate an index- X DAE on Ai or an index- K + 1 DAE on a symplectic embedding 
of Ai so as to preserve the constraints and i_m k uj. 

Finally, we mention another class of integrators for the holonomic case, known as 
spark, for Symplectic Partitioned Additive Runge-Kutta [9]. These generalise rattle 
to higher order. They are partitioned (use different RK coefficients for the q and p 
components) and additive (use different RK coefficients for the constraint and regular 
forces). The holonomic assumption is used in two critical steps: first, it means that 
the flow of the constraint force is given by Euler's method; second, it means that the 
g-component of the constraint forces vanishes. This allows their RK coefficients to be 
arbitrary, which means that the RK coefficients of the p-component can be arbitrary, 
and can be chosen to include stages at the endpoints. Thus, this approach does not 
immediately give an algorithm for problems of the type (V,Ai',H) mentioned above. 
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The situation is similar to the relationship between splitting methods and RK methods; 
we do not know if spark can be adapted to more general constraints. 
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